Forces between elongated particles in a nematic colloid 
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Using molecular dynamics (MD) simulations we study the interactions between elongated colloidal 
particles (length to breath ratio ^ 1 ) in a nematic host. The simulation results are compared to 
the results of a Landau-de Gennes elastic free energy. We find that depletion forces dominate 
for the sizes of the colloidal particles studied. The tangential component of the force, however, 
allows us to resolve the elastic contribution to the total interaction. We find that this contribution 
differs from the quadrupolar interaction predicted at large separations. The difference is due to the 
presence of nonlinear effects, namely the change in the positions and structure of the defects and 
their annihilation at small separations. 
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I. INTRODUCTION 

Liquid crystal colloids belong to a special class of col- 
loidal systems. Long-range orientational order of the liq- 
uid crystal molecules gives rise to additional long-range 
interactions between the colloidal particles The 
presence of defects, due to topological restrictions, con- 
trols the symmetry of this interaction, that may be of 
dipolar or quadrupolar symmetry H, Clustering, 
superstructures, and new phases are immediate conse- 
quences of additional anisotropic interactions 0, Q, I3| ■ 

The long-range forces between particles of any shape 
can, in principle, be calculated using direct integration 
over the director field 0- Ready-to- use expressions 
are available for spherical particles |lfll | , two-dimensional 
disks ^J, etc. For smaller particle-particle separations 
nonlinear effects from the elastic free energy come into 
play. Minimization of the Landau-de Gennes free energy 
with respect to the tensor order parameter can be used 
to take into account the relative position of the defects 
and the variation of the nematic order parameter around 
the colloids For even smaller separations, deple- 

tion forces, density variation, and presmectic ordering 
of a nematic liquid crystal near colloidal particles can- 
not be ignored. Then the density functional approach 
|l2l |l 4l[ or^ alternatively, computer simulation tech- 
niques |15lll6l[l7l| can be used. 

In this paper we study the interaction between elon- 
gated colloidal particles suspended in a nematic liquid 
crystal. The liquid crystal molecules are taken to be 
homeotropically anchored, that is, their preferred orienta- 
tion is normal to the colloid surface. We compare the re- 
sults of two methods: molecular dynamics simulation and 
minimization of the phenomenological Landau-de Gennes 
free energy. 

It has already been shown theoretically and using 
MD simulation that isolated small elongated par- 
ticles, with homeotropic boundary conditions, minimize 



the free energy by orienting perpendicularly to the di- 
rector. Thus, we consider the particle symmetry axes 
parallel to each other and perpendicular to the far field 
director. In this geometry, the defect structure around a 
single particle has been studied in detail |3| • The config- 
uration with two 1/2 strength defects is preferable ener- 
getically, giving rise to a quadrupolar-type of interaction 
between two such particles |0. Analytical expressions 
for the long-range forces derived from the Frank elastic 
free energy exist, as well as numerical studies based on 
the tensor order parameter formalism 0, |23| • 

However, more detailed studies of the forces showed 
that the interaction between the particles is no longer 
quadrupolar at small separations [ij , due to a change in 
the relative position of the defects. It was shown that the 
long-range repulsive interactions can become attractive 
for small anchoring strengths while remaining repulsive 
for all orientations for strong anchoring. As the distance 
between the particles decreases, their preferred relative 
orientation with respect to the far field nematic direc- 
tor changes from oblique (7r/4 for the pure quadrupolar 
interaction) to perpendicular. 

In this paper we use MD simulations to confirm this 
conclusion at even smaller separations, where the deple- 
tion forces play an important role. To resolve the elas- 
tic contribution to the total force we measure separately 
the normal and the tangential component of the force. 
The normal component is much larger than the tangen- 
tial component and it is practically unaffected by the 
relative position of the defects. The tangential compo- 
nent, however, has a dependence on the particle separa- 
tion that is qualitatively the same as that predicted by 
the minimization of the Landau-de Gennes free energy. 

The paper is organized as follows. In Sec. ^ we de- 
scribe the geometry, the molecular model used for MD 
simulation, and the technique of minimization with adap- 
tive meshing. The director orientation, order parameter, 
and density maps, as well as forces between the parti- 
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FIG. 1: Geometry (yz cross-section is shown): two cylinders 
of radius R are immersed in a liquid crystal host. The host is 
modeled as a solution of Gay-Berne particles and the orienta- 
tion of the director is fixed at the top and bottom walls of the 
simulation box [z = ±Lz/2). Periodic boundary conditions 
are applied in the x, y directions. 




cles are presented in Sec. IIIII In Sec. IIVI we make some 
concluding remarks. 

II. MOLECULAR MODEL AND SIMULATION 
METHODS 

The geometry considered in this work is shown in Fig- 
ure n Two cylindrically shaped colloidal particles of ra- 
dius R are immersed in a liquid crystal. The particles are 
separated by a distance d measured from their symmetry 
axes. The director orientation is fixed at the top and the 
bottom walls, parallel to the z axis and perpendicular to 
the symmetry axes of the cylinders. Boundary conditions 
ensure that the director n far from the colloidal particles 
is parallel to the z axis. The rod length is considered to 
be infinite: in the simulations this simply means that the 
cylinder spans the x-dimension of the periodic box. 

A. MD simulation 

Molecular dynamics simulations were carried out using 
a soft repulsive potential, describing (a ppr oximately) el- 
lipsoidal molecules of elongation k = 3 [T3 . The systems 
consisted of iV = 8,000 and 64,000 particles (see table 13 
for details). A reduced temperature fceT/eo — 1 was 
used throughout. The system size was chosen so that 
the number density is pa^ = 0.32, within the nematic 
phase for this system. Here ctq is a size parameter, eg is 
an energy parameter (both taken to be unity). 

The system is confined in the z direction, to provide 
uniform orientation of the director far from the parti- 
cles along the z axis. The interaction of molecule i with 
the colloidal particle (rod) and the wall was given by 
a shifted Lennard- Jones repulsion potential having ex- 
actly the same form as in Refs. Iialigj. This provides 
homeotropic orientation with a strong anchoring of the 
molecules at the wall. 



FIG. 2: MD simulation results: snapshot of the system. Num- 
ber of particles A'^ = 8000, colloid radius R/ao — 3, colloid 
separation d/ao = 10, a = 0. Color coding emphasizes parti- 
cle orientations. 



The radius and the length of the rod were steadily 
increased from zero to the desired value during 10"^ steps. 
Then the system was equilibrated for 10^ steps. During 
equilibration we scaled the velocities of the molecules to 
achieve ksT/eo = 1. An equilibrated snapshot of the 
system is shown in Fig. [3 

The production run for every angle a was 10^ steps. 
The force F on the rod was calculated using the repulsive 
force fi from the rod on the particle i 

N 
i=l 

The local tensor order parameter Q{r) was calculated 

as 

1 (3 11 

Qap{zi,yj) = ^J— ^ < - (UkaUkp) - -^Sap > , (2) 

where there are n^i jy molecules present in each bin {i,j}, 
da[3 is the Kronecker delta, (■••) denotes an ensemble 
average, a,(3 = x,y,z. Diagonalizing the Qaf3 tensor, 
for each bin, gives three eigenvalues, Qi, Q2, and Qa, 
plus the three corresponding eigenvectors. The eigen- 
value with the largest absolute value defines the order 
parameter S for each bin. The biaxiality B is calculated 
as the absolute value of the difference between the re- 
maining two eigenvalues of the tensor order parameter. 
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TABLE I: Systems studied 



Colloid radius 


Particles 


Box size 


Production 




N 


x/cro, y/(JQ, z/ao 


(steps, 10^) 


3 


8,000 


10, 50, 50 


10 


5 


64,000 


20, 100, 100 


5 



B. Landau-de Gennes free energy minimization 

In the framework of the continuum theory, the sys- 
tem can be described by the Landau-de Gennes free en- 
ergy [m 



T{Q} 



HQ) 



(3) 



where / [Q) is a function of the invariants of Q, the 
symmetric tensor order parameter, and the integration 
extends over the sample volume. Here we adopted the 
one-constant approximation for the elastic free energy 
(decoupling of the spatial and spin rotations). 

To simplify the calculations, we used a two- 
dimensional representation of the tensor order parame- 
ter. This representation is justified for a uniaxial nematic 
with the director in the yz plane. Indeed, as we will see, 
the MD simulation results show that the director is con- 
strained in the yz plane. The nematic is uniaxial in the 
bulk and is slightly biaxial in the defect core. 

In this approach it is possible to rewrite Q in terms 
of components of the director n and the scalar order pa- 
rameter Q 



Qij = Q niUj - -Stj 



(4) 



Then the Landau-de Gennes free energy density reads 



fiQ) 



[Trg2 



(5) 



where a is assumed to depend linearly on the tempera- 
ture, whereas the positive constant c is considered tem- 
perature independent. Note, that the invariant that cor- 
responds to the cubic term of the tensor order parameter 
vanishes in the two-dimensional nematic. 

For this free energy the liquid crystal in the nematic 
state has the order parameter Q^q = \j2ajc. The elastic 
constant L is related to the Frank elastic constant by 
K = ALa/c. In our calculations we used a — l,c = 2 
which gives Qeq = 1. 

The free energy J-{Q} was minimized numerically, us- 
ing finite elements with adaptive meshes. The geometry 
used in the numerical calculation is shown schematically 
in figure n During minimization we used a square inte- 
gration region fl of size 40i? x 40i?, where R is the radius 
of the colloidal particle. The area 17 was triangulated us- 
ing a BL2D subroutine i23 | . The tensor order parameter 
Q was set at all vertices of the mesh and was linearly 



interpolated within each triangle. Using standard nu- 
merical procedures the free energy was minimized under 
the constraints imposed by the boundary conditions, i.e., 
strong homeotropic anchoring at the particle perimeters 
and uniform alignment at the outer boundary. 

Finally, a new adapted mesh was generated iteratively 
from the result of the previous minimization. The new 
local triangle size was related to the free energy variations 
of the previous solution, assuring a constant numerical 
weight for each minimization variable. The final meshes 
with a minimal length of ^ 10~^R had about 2 x 10'' 
minimization variables. 

The tangential component of the interparticle force, 
Ft , was calculated numerically from the free energy, = 
—r^^dT /da. 



III. RESULTS 



A. Defect structures 



A typical director orientation together with the den- 
sity, order parameter, and biaxiality maps are shown in 
Fig. 13 For large separations, a pair of defect lines forms 
perpendicular to the director, next to each particle, and 
stays at the same positions as in the case with a single 
colloidal particle. The director distortion vanishes very 
quickly in the liquid crystal bulk, and the core region 
extends over a few molecular lengths [l^ . 

On reducing the separation, at a certain value the posi- 
tions of the defects start to change, as well as their inner 
structure (see Fig.|3|D,c). Finally, when the particles are 
about to merge, two out of four defects vanish to ensure 
that the total topological charge in the system is zero 

(Fig. EH). 

This scenario agrees qualitatively with the results 
of the minimization of the Landau-de Gennes free en- 
ergy Using MD simulation, we are also able to 
observe the annihilation of the defects, which is not ac- 
cessible on length scales where the phenomenological ap- 
proach is applicable. From the density map shown in 
Figure EJi one can already foresee the strong influence of 
depletion effects at small separations. The order parame- 
ter map (Fig. (Jb) illustrates the dislocation of the defects 
and the changes in the defect cores. 
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FIG. 4: MD simulation results: components of the force paral- 
lel to the particle-particle separation vector, Fy , as a function 
of the particle-particle separation d. Colloid radius R/ao = 3. 



FIG. 3; MD simulation results. Maps of (a) density at the 
separation d/ao = 16, (b) order parameter at the separation 
d/ao — 13 , and (c) biaxiality at the separation d/ao ~ 13. 
(d) shows the map of the Uy component of the director at 
the separation d/ao = 6. The yz cross section of the system 
is shown. Rod radius R = Sao- The director far from the 
particle is constrained along the z axis. 



B. Force between the particles 

Figure 01 shows typical force curves as a function of 
the particle-particle separation d for colloidal particles of 
size R/ctq = 3. From the parallel component of the force 
(projection of the force on the particle-particle separa- 
tion vector) one can see that the depletion forces indeed 
dominate; the force has oscillations due to the density 
modulation of the liquid crystal close to the particle sur- 
faces (presmectic ordering). More detailed analysis of 
the force curves, shown in the inset of Fig. 01 reveals 
that there is a small difference in the decaying tails of 
the force curves, which can be attributed to the elastic 
contribution of the order parameter field/defects around 
the colloidal particles to the total interparticle interac- 
tion. 

To emphasize the contribution of the elastic force, we 
plot the tangential component of the force (perpendicu- 
lar to the particle-particle separation vector) in Fig. \^ 
We first look at the situation with a = tt/A. It is clear 
that there is a non-zero tangential component, which de- 
cays with the distance, i.e. there is a non-zero force which 
tends to align the particle-particle separation vector per- 
pendicular to the director (or at some angle which is less 
than 7r/4). This is already different from the prediction 
of the quadrupolar interaction, where the minimum of 
the free energy is at a = 7r/4, and agrees with the re- 
sults of the Landau-de Gennes free energy minimization 



for small particle separations, see Ref. [llj. 

The situation with a — 0, 7r/2 is somewhat different: 
there is a scatter in the value of the tangential component 
of the force, accompanied with a large error in the mea- 
surements. Analysis of the configurations suggests that 
this is due to the degeneracy in the possible defect posi- 
tions. Indeed, if the particles move close to each other, 
the defects change their positions. The vector between 
the defects belonging to the same particle tilts with re- 
spect to the director. If a = or tt/2, the tilt angle 
can be either positive or negative: both configurations 
are equivalent and have the same energy. However, there 
is a barrier between these two configurations. For small 
particles, this barrier is of the order of k^T, and the 
defects can switch between two equivalent configurations 
during the simulation run. This 'drift' of the defects does 
not affect the parallel component of the force, but leads 
to the 'averaging' of the tangential component to zero. 
In addition, this drift is rather slow on the timescale of a 
molecular simulation, and provides a large scatter in the 
value of the tangential force. 

To overcome this problem, we also studied particles 
with a bigger radius, R — Scto- The elastic contribution 
to the force is larger in this case. As a consequence, the 
energy barrier between the two equivalent configurations 
is also larger. It is harder for the defects to overcome 
this barrier, which becomes larger than ksT for this size 
of the colloids. This is clearly seen from the tangential 
component of the force plotted in Fig.|Hl the scatter of the 
data is much smaller and there is a clear decrease with the 
decrease of the particle-particle angle. One can also see 
that there is a non-zero tangential force for a = tt/2, 7r/4, 
i.e. the particles tend to align in a way that their center- 
center separation vector is perpendicular to the director. 

Finally, we show the results of the Landau-de Gennes 
theory in Fig. The tangential component of the force 
is plotted as a function of the interparticle distance d, for 
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FIG. 5: MD simulation results: component of the force per- 
pendicular to the particle-particle separation vector, Ft, as a 
function of the particle-particle separation d. Colloid radius 
R/ao = 3. Error bars for a — n/2 are too big to show on the 
plot. Smooth curves are to guide the eye. 



Gennes free energy, to study the interaction of two elon- 
gated colloidal particles embedded in a nematic host. 
Our results show that the particle-particle interaction is 
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FIG. 6: MD simulation results: component of the force per- 
pendicular to the particle-particle separation vector. Ft, as a 
function of the particle-particle separation d. Colloid radius 
R/ao = 5. Smooth curves are to guide the eye. 



FIG. 7: Landau-de Gennes theory: component of the force 
perpendicular to the vector joining the centers of the colloids 
as a function of the colloid separation, d, at three different 
orientations {a = 0, 7r/4, 7r/2). d is in units of the nematic 
correlation length ( = \fTJa. a) (,/R^ 0.2; h) C/R^ 0.633 



no longer quadrupolar at short distances due to a change 
in the relative position of the defects. MD simulation 
results also show that for small particles the depletion 
force dominates but contributes mostly to the interparti- 
cle radial force. The tangential contribution to the force 
is of elastic origin. Its dependence on the center-center 
separation is in qualitative agreement with the results of 
the free energy minimization. 



three different orientations of the particle-particle sepa- 
ration vector, a = 0,7r/4, 7r/2. 

Comparing the results of the Landau-de Gennes the- 
ory, Fig.[7| to the MD simulation results, Fig.|Hl o ne se es 
qualitative agreement for C,/R = 0.2 where Q — \/ L/a is 
the nematic correlation length. 

IV. CONCLUSIONS 

We used two independent techniques, molecular dy- 
namics simulation and minimization of the Landau-de 
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